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ABSTRACT 


A local geoid model to predict the geoid heights in the vicinity of Monterey 
Bay, California, was developed to use Global Positioning System (GPS) differential 
positions and known Mean Sea Level (MSL) with the method of collocation. The 
local geoid models were based on Rapp's 360 degree x 360 order global geoid 
model determined from gravity measurements. Control data were adjusted by least 
squares to solve for the parameters in the local geoid model. Also studied were 
factors that affected the GPS-measured ellipsoid height differences. These included 
(1) comparing GPS differencing solutions, (2) standard error of GPS observations, 
(3) corrections for surface meteorological values, and (4) observation durations for 
GPS. 

The data used in this research were taken from GPS measurements on the 
campus of the Naval Postgraduate School (NPS), an area about 100 m x 630 m and 
in an area approximately 15 km x 33 km near Monterey, California. The time 
period was from February 5, 1988, to May 12, 1988. 

The accuracy of the predicted geoid heights is + 2 cm if a six-parameter model 
is used for the large area, and + 2 to 10 mmif a five-parameter model is used for the 


NPS campus. 
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I. INTRODUCTION 


The Global Positioning System (GPS) is able to establish precise relative 
positions in the World Geodetic System of 1984 (WGS 84). A Trimble GPS 
receiver, which has the capability of measuring carrier phase, was used to 
determine the vector base line in space. The components of the base line are 
expressed in terms of cartesian coordinate differences (AX, AY, AZ) [Remondi, 
1984]. These vector base lines can be converted to distances, azimuths and the 
ellipsoid height differences, Ah, relative to the WGS 84 Ellipsoid. 

The results of several tests and operations have clearly shown that GPS 
Survey methods can replace conventional horizontal survey methods. 
Comparable accuracies have also been achieved for GPS-derived ellipsoid height 
differences, Ah. The problem of converting these ellipsoid height differences, 
Ah, to orthometric height differences, AH, remains to be resolved. For example, 
in engineering surveying applications WGS 84 coordinates must be transformed 
to the North American Datum of 1983 (NAD 83) system. The GPS obtains 
ellipsoid heights, rather than the orthometric heights; the geoid height, N, must 
be calculated to obtain the latter. One of the problems in this transformation 15 
the accurate determination of the local geoid height, N. 

For GPS survey applications, a geoid model should provide geoid heights 
with an accuracy commensurate with that of the ellipsoid height, h, so the 
accuracy of the derived orthometric heights is not reduced. In the future 
differential positioning will be widely used in GPS surveying, and, therefore, 
only the geoid height differences between stations will be required. 

Geoid height computation techniques include the following: 

(1) Geoid height differences in the U.S. can be determined from gravity 
data and the Stokes’ integral method, or from astrogravimetric data and least 
squares collocation methods. These methods lead to uncertainties that are 
typically 1 to 10 cm for distances of less than 20 km and 5 to 20 cm for distances 
between 20 to 50 km [Zilkoski, 1988] 

(2) A comparison of a data set from GPS with gravimetrically determined 
geoid heights using least squares collocation techniques shows discrepancies 
between the two data sets of about + 2 cm for a maximum intersection distance of 
approximately 50 Кт [Denker and Wenzel, 1987]. 


(3) Mean gravity anomalies, deflections of the vertical and a geopotential 
model calculated to degree and order 180 have been used to determine geoid 
heights in the area bounded by (34? € 9 € 42?, 18? € À € 28?) [Tziavos, 1987]. 
The method used was that of least squares collocation. By using empirical 
covariance functions for the data, suitable weighting functions for the different 
sources of observations, and the optimum cap radius around each point of 
elevation, an accuracy better than + 0.60 m was obtained for geoid heights. 

The main objective of this thesis is to use a model that predicts N in a region 
near Monterey, California, from the GPS differential positions and known mean 
sea level (MSL) using the method of collocation. Also studied were factors that 
affected the ellipsoid height differences, Ah, obtained from GPS measurements, 
such as comparing GPS differencing solutions, standard error of GPS 
observations, corrections for surface meteorological value , and observation 
durations for GPS measurements. The data used in this research were taken 
from GPS measurements on the campus of the Naval Postgraduate School (NPS) 
an area about 100 m x 630 m and in an area approximately 15 km x 33 km near 
Monterey, California. The time period was from February 5, 1988, to May 12, 
1988. 

The accuracy of the predicted geoid height is + 2 cm if a six-parameter 
model is used for the large area, and + 2 to 10 mm if a five-parameter model is 
used for the NPS campus. 


П. GEOID HEIGHT 


A. THE GEOID 

Surveyors and engineers, in most cases, are interested in the orthometric 
height, H, as measured above the reference surface of the geoid. The ocean is 
considered to be freely moving, homogeneous and only subject to the force of 
gravity. When a state of equilibrium is achieved, the surface of this idealized 
ocean assumes a level surface of the gravity field. It may be regarded as also 
extending under the continents. This level surface is called the geoid. If the 
potential, W, is given as a function of the coordinates r, ф and A, then the geoid is 
given by [Moritz, 1984] 


W - W(r$,À) = №, 


The geoid is a closed and continuous level surface which extends partially 
inside the solid body of the earth. The direction of the gravity vector at any 
point (plumb line or vertical) is normal to the geoid. The curvature of the geoid 
displays discontinuities at abrupt density variations. Consequently, the geoid is 
not an analytic surface, and therefore not a practical reference surface for 
position determinations. The geoid however, is well suited as a reference surface 
for potential or height differences, which are obtained by precise levelling in 
combination with gravity measurements. 

To establish the geoid as a reference surface for heights, the ocean water 
level is recorded and averaged over long intervals (2 1 year) using tide gauges. 
The MSL thus obtained represents an approximation to the geoid. The National 
Geodetic Vertical Datum of 1929 (NGVD 29) was derived for land surveys from 
a general adjustment of the first order levelling net of both United States of 
America and Canada. In the adjustment MSL was observed at twenty-one tide 
stations in the United States and five in Canada. The geoid established by this 
method may deviate by + 1 to + 2 m from a level surface due to periodic, 
nonperiodic and secular variations [Torge, 1980]. 


B. THE WGS 84 ELLIPSOID 

The development of WGS 84 [DMA, 1987] was initiated by the United States 
Department of Defense for navigation and weapon systems. The Defense 
Mapping Agency (DMA) developed WGS 84 as a replacement for WGS 72. The 
defining parameters and reference frame orientation of the WGS 84 Ellipsoid 
and the WGS 84 Ellipsoid Gravity Formula are those of the internationally 
sanctioned Geodetic Reference System of 1980 (GRS 80). Accordingly, a 
geocentric equipotential ellipsoid is defined by the semimajor axis (a), the 
flattening (f), the equatorial gravity (ya), and the angular velocity (œ). The WGS 
84 Ellipsoid used the values: 


= 6378137 т 

f =  1/298.257223563 

Ya =  987.03267714 gals 
© =  7292115x10? rad/s 


The reference system for GPS is WGS 84. The precise geocentric 
coordinates obtained from GPS receivers are in WGS 84. 


C. THE ORTHOMETRIC HEIGHT, H, THE ELLIPSOID 


HEIGHT, h AND THE GEOID HEIGHT, N 

In geodetic applications three different surfaces or earth figures are 
normally involved. First is the earth's actual topography; second, the geometric 
surface, or ellipsoid and; third, the equipotential surface, the geoid. The 
relationship between the earth's topography, the ellipsoid and the geoid in a 
section through the earth's surface is shown in Figure 1. 


NORMAL TO ELLIPSOID 


NORMAL TO GEOID 





Figure 1. Relationship between the earth's surface, the geoid and 
the ellipsoid. 


Two features in the figure are of particular interest : 

1. The deflection of the vertical, 0, defined by Pizzetti as the angle at the 
geoid between the direction of the plumb line (normal to geoid) and 
the normal to the ellipsoid through the point, P, on the geoid [Torge, 
1980]. | 

2. The vertical separation, N, between the geoid and the ellipsoid. 


The deflections of the vertical and the geoid heights, N, depend on the 
ellipsoidal coordinates and, hence, on the parameters of the reference ellipsoid 
and its position with respect to the earth. If they are referred to the 
geocentrically situated mean earth ellipsoid, then they are referred to as absolute 
quantities; otherwise, they are relative quantities. The absolute deflections of the 
vertical in flat terrain and the highlands assume values between one and ten 
seconds of arc; in mountainous areas, they vary between one-half and one minute 
of arc. As a result of density variations within the earth, the geopotential 
surfaces, including the geoid, have irregular shapes. Absolute geoid heights, 
however, rarely exceed 100 m [Torge, 1980]. 

The orthometric heights, H, shown in Figure 2 are referred to an 
equipotential surface, the geoid. The orthometric height of a point on the earth's 
surface is the distance from the reference surface to the point, measured along 
the plumb line normal to the geoid. The ellipsoid height, h, of a point is the 
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distance from the reference ellipsoid to the point, measured along the line which 
is normal to the ellipsoid. For purposes of simplicity, the orthometric height, H, 
and the ellipsoid height, h, are shown to be along a common vertical. In most 
cases, this would cause a very small error that is considered insignificant 
compared to present uncertainties of the geoid height N estimates. The geoid 
height, N, is defined: 


vn 
: 


“sa TEA AA 


ELLIPSOID 





Figure 2. Relationship between the geoid height, N, the ellipsoid 
height, h and the orthometric height, H. 


The orthometric height has greater physical meaning than the geometrical 
ellipsoid height. The orthometric height has traditionally been determined by the 
technique of levelling in which increments of height are obtained from the 
intersection of the line of sight of a level instrument, tangential to the 
geopotential surface and passing through the level axis, with two graduated rods. 
Accurate orthometric height information is needed for precise engineering 
operations such as the construction of dams, pipelines, and tunnels. 

Geoid heights have been accurately determined for some major geodetic 
datums such as the North American, European and Australian. These datums are 


well supplied with astrogeodetic deflections and have fair gravity coverage. The 
standard error of relative geoid heights in these areas is about 2 or 3 m. 


D. THE GEOID HEIGHT FROM GRAVITY MEASUREMENT 

Both the U.S. National Geodetic Survey (NGS) of National Oceanic and 
Atmospheric Administration Service (NOAA) and Trimble Navigation versions 
of Rapp's 360 degree x 360 order model (OSU86F) were available to me. Rapp 
computed two potential coefficient fields that are complete to degree and order 
360 [Rapp, 1986]. One field (OSU86E) excludes geophysically predicted 
anomalies, while Rapp's other model (OSU86F) includes such anomalies. These 
fields were computed using a set of 30-minute mean gravity anomalies derived 
from satellite altimetry in the ocean areas and on land from standard 
measurements. 

Gravity anomalies can be observed and then used to compute the geometric 
deviation of the geoid from the ellipsoid. The expression for the gravitational 
potential is written in the following form [Rapp, 1986]: 


re) = eM 1 + У Bi y ( €, cosm +5 sinmA | ЖС | 


Y m=0 


where 
т, 6, A : geocentric coordinate 
kM : geocentric gravitational constant 
a  : equatorial radius of the reference ellipsoid 


C_»5,,,: fully normalized potential coefficients 
P „ : fully normalized Legendre function of degree | and m 


Y  :normal gravity 
The potential at a point, U, is the scalar sum of V and centrifugal force potential 
V' [Ewing, 1976]: 


О(т,ф,А) 2 V(r,$,A) + V'(r,6,A) 


The difference between observed gravity potential W and the computed normal 
gravity potential U is denoted by T [Moritz, 1984], so that 
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W(r,$,À) = UG,o,A) + Т(т,ф,А) 

compared the geoid defined by the potential W3 is given by 
W(r,0,A) = Mic 

A reference ellipsoid with the same potential, W, = U „is given by 


О(т,ф,А) = W, 


A point P on the geoid is projected onto the point Q on the ellipsoid along the 
normal to the ellipsoid PQ = N. PQ is the distance between geoid and ellipsoid at 
the point and is called the geoid height (Figure 3). 


ELLIPSOID 
U = Wo 





Figure 3. Geoid and ellipsoid (Moritz [1984, Fig. 2-12]). 
The disturbing potential T can be written as [Rapp, 1986]: 
КМ > ] l A — 
TA) = T^y [f| X € C YR 
127 m=0 a=0 
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where 


Ce = (C_,a=0, and 5 ,а=1) 


Y (6,2) = (P, (sing)cosma, a=0, and Р (ѕіпф)ѕіптл, ade 


Since the N is relatively small compared to the ellipsoid geocentric radius R', 
then [Ewing, 1976] 


and 


where y is the normal gravity at (@, A). 
The gravity anomaly, Ag, in the Molodensky surface free-air anomaly sense 
is given by [Moritz, 1984]: 


Agp = Ep - Yo 


The boundary condition that relates Ag and T is [Moritz, 1984] 


where h is the distance along the plumb line direction. Neglecting deflections of 
the vertical we have [Rapp, 1986] 


k = е 1 -а A x V 
п ООО мю мш шуу) 
Y 1-2 0 


Г т= 0 a= 





Ag = 


where e (i= h, y) are ellipsoidal corrections. 
The Stokes function, S(y), can then be used to solve the geoid heights above the 
geodetic ellipsoid [Ewing, 1976]. 


1 21 іл : 
N, = ae f Ag(y,o)S(y)sinydwydo: 


where 
Ү, “the mean value of normal gravity 


: the angular distance between the point where N is being 
determined and the area where the effect of Ag is being 


considered 
a: the azimuth from the affected point to that causing the effect 


Ag  :the gravity anomaly 
S(y) : the Stokes function 


The Stokes function is given by [Ewing, 1976] 


S(y) = ese T Lom Zu 5 cosy - 3 cosy In(sin 2-4 sin” 2) 


The gravitational potential using degree 360 and order 360 has been tested 
through comparison of Doppler station geoid heights with geoid heights from 
Rapp's versions of geopotential models. The agreement between the two geoid 


height measurements is approximately + 1.6 m [Rapp, 1986]. 
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ПІ. GEOID HEIGHT FROM GPS 


A. THE GLOBAL POSITIONING SYSTEM 

The Global Positioning System (GPS) calls for a precise navigation system 
divided into three segments: space segment, control segment and user equipment. - 
The space segment will consist of three orbital planes of satellites at inclinations 
of 120° in circular orbits at altitudes of 20,000 km. Each plane will eventually 
contain six to eight satellites to give the three dimensions of position, velocity and 
precise time 24 hours a day anywhere in the world. The control segment consists 
of the ground stations necessary to track the satellites, monitor the system 
operation, and periodically provide corrections to the navigation and time 
signals. Each satellite broadcasts signals containing information on its position. 
The GPS satellite transmits signals at two L-band frequencies (1227 and 1575 
MHZ) to permit corrections for ionospheric corrections. The signals are 
modulated with two codes: P, which provides for precise measurement, and C/A, 
which permits easy lock-on to the desired signal. The user segment consists of 
the equipment necessary to convert the satellite navigation message into useful 
navigation information. 

The navigation message contains the data that the user's receiver requires to 
perform the operations and computations for successful navigation with GPS. 
The data includes: (1) information on the status of the Space Vehicle (SV); (2) 
time synchronization information for the conversion of the C/A to P code; (3) 
parameters for computing the clock correction and the ephemeris of the SV; and 
(4) corrections for delays in the propagation of the signal through the 
atmosphere. In addition the data contain almanac information to define the 
approximate ephemeris and to give the status of all the other SV information 
which is required for use in signal acquisition [Milliken, 1980]. Ranges to the 
satellites are determined by signal transit times multiplied by the speed of light 
(299,792,458 m/sec). The transmitted message contains ephemeris parameters 
that enable the user to calculate the position of each satellite at the time of the 
transmission of the signal. 
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B. CARRIER PHASE MEASUREMENTS 

GPS measurements can be made using the pseudo-range and the carrier 
phase. The pseudo-range is essentially a measurement of distance contaminated 
by clock error. When four satellites are observed simultaneously, the three 
dimensional position of the ground receiver can be determined along with the 
receiver clock offset at a single epoch. The accuracy of pseudo-ranges is affected 
by multipath effects, which depend on the antenna design, and its height above the 
ground. 

Carrier phase measurements are more precise than pseudo-range and are 
not as vulnerable to multipath effects. They can be used to compute the precise 
base line components AX, AY, AZ between two receivers. Phase measurements 
are made by beating the received carrier with the signal from a local oscillator 
internal to the GPS receiver. The slant range from a GPS receiver to a satellite 
can be modelled in terms of time. It takes the signal time to travel between the 
satellite and the receiver or the equivalent number of cycles. The range of the 
cycles will consist of an integer and fractional number of cycles. When a 
receiver locks onto the carrier signal, it can immediately measure the fractional 
part and begin counting subsequent integer cycles, but it can not measure or 
account for the initial integer number of cycles that preceded the initial fractional 
part. The initial integer ambiguity which biases the subsequent measurements is 
called the initial integer ambiguity bias. 

GPS uses a one-way carrier beat phase. The GPS satellite and receiver are 
controlled by separate clocks. The satellite clock generates the signal, and the 
receiver clock detects when the signal arrives. An error in the synchronization 
of the clocks of 1 microsecond creates an error in range of 300 m. 


C. ONE-WAY CARRIER PHASE MEASUREMENT 
DIFFERENCING 


A single difference, SD(j,i), is formed by differencing carrier beat phase 
observables from two receivers 1 and 2 at the same observation epochs i of same 
satellite j}. The equation is given by [Remondi, 1984]: 


SD(ji) 2» S(Q,j) - S(1,,1) 
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where S(k,j,i) is the raw, unpreprocessed, fractional phase plus the count made at 
epoch i by receiver k for satellite j}. The main advantage of the single difference 
is that it reduces or eliminates satellite orbital and clock errors, because they are 
common to both receivers. Its disadvantage is that one can not exploit the integer 
nature of integer ambiguities. Thus, for short base lines, the ultimate in accuracy 
may not be achievable [Remondi,1985]. 

A double difference, DD(j,k,i), is formed by differencing single differences 
between a reference satellite j and another satellite k at the same epoch ı. The 
equation is given by [Remondi, 1984]: 


DD(j,k,i) = SD(k,i) - SD(j,i) 


The advantage of the double differences is that the receiver clock dependent 
terms are eliminated because the differences for each epoch are correlated. The 
significance of the removal of those terms is to reduce from nanoseconds to 
microseconds the timing accuracy required to achieve one cycle accuracy. For 
short base lines, the integer ambiguities can be isolated, since the contribution 
made by the clock drift is reduced [Remondi, 1985]. The Trimble 4000SX 
receiver achieves sub-microsecond accuracy by using the C/A code timing 
information [Ashjaee, 1985]. 

A triple difference, TD(j,k,i), is formed by differencing the double 
differences for the same satellite pair at some integer of succeeding epochs i+] 
[Remondi, 1984]. 


TDg,k,1) = DDG,k,i+1) - DDG,k,i) 


The advantage of the triple difference is that it eliminates all the time independent 
terms, namely the initial integer ambiguities, and becomes insensitive to the 
initial ambiguities and any cycle slips when.the receiver loses lock. The 
disadvantage of the triple difference is, another level of correlation, loss of 
resolution and a greater number of observations. Triple differences are already 
correlated with respect to satellite because of the underlying double differences 
and are further correlated with respect to time because consecutive triple 
observations will have common DD(j,k,i+1) terms. 
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For short base lines local area integer ambiguities can easily be resolved 
because unmodelled errors are highly correlated between the two antenna sites and 
are mostly eliminated by differencing. Algorithms can take advantage of the 
integer nature of the initial ambiguities and solve for them [Remondi, 1984]. 


D. GEOID HEIGHT FROM GPS 

Surface fitting techniques can be used with GPS-derived geoid heights. GPS 
Stations are likely to be close together, of the order of a few tens of km, and the local 
geoid can be estimated directly if levelling data is available in the area. This 
together with it's relative simplicity makes the method practical for correcting GPS 
heights. 

It is possible to use the two sets of elevations (that is, levelled and GPS- 
derived) to define two distinct planes. The published levelled elevations are 
referred to the geoid and the GPS elevations are referred to an ellipsoidal surface. 
If both elevations are made equal at one bench mark (control station), then in 
general, the other bench marks will have two elevation values. After several 
different models were studied two mathematical surface models were chosen in this 
study. The five-parameter model found suitable for use on the NPS campus is 


No(X,Y,Z) = hy + Ah-H + aAY + bAX? + cAY? + dAXAY 


and the six-parameter model selected for a larger area near Monterey, California, is 
given by 


N,(X,Y,Z) = ho + Ah- H + AX +D'AZ + CAX? + d'AY? + e'AXAY 


where 

N,(X, Y,Z) : the global geoid height, obtained by gravity measurements 

họ : the ellipsoid height of the reference point, including a constant 
correction to N,(X,Y,Z) at the reference station 

Ah :the ellipsoid height difference with respect to the reference station 

H :the published or levelled elevations at each control station 

a, b, c, d, a', b', c', d', e' : the coefficients to be determined 

AX, AY, AZ : the coordinate differences in WGS 84 
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These models can be solved to determine the east-west and north-south tilts 
which are absorbed by two rotations, one around the north axis (AY) and the 
other around the east axis (AX) in the horizon system. Their separation is 
absorbed by the scale correction [Zilkoski, 1988]. The local geoid heights, N, 
can be found from the equation 


N = N,(X,Y,Z) + AN 

where AN is the variation of geoid height in local area. The equation is given by 
AN = H+N,(X,Y,Z) - Ah 
AN = + аду + bAX? + cAY? + dAXAY 


or 
AN 


hy + a'AX + PAZ + c'AX? + d'AY? * e'AXAY 


To solve for hg, а, Б, с, 4, а, Б, с, 4 апд е, the global geoid heights 
N,(X,Y,Z) can be obtained from the Geoid.exe program described in Chapter 


IV. The geoid model can be rearranged to give 


H = h, + Ah -N,(X,Y,Z) + aAY + bAX? + cAY? + dAXAY 


OT 
H = h, + Ah - N (X,Y,Z) + a'AX * D'AZ + c'AX? + d'AY? + e'AXAY 


Then ho. a, b, c, d, a, b', c', d' and e' can be solved by the least squares method. 


The geoid height, N, found by this method appears to be adequate for areas up to 
50 km x 50 km where the geoid is smooth [King, 1985]. 


E. METHOD OF COLLOCATION 

The method of collocation was derived from least squares interpolation. 
This method was used to predict the geoid height, N, for a local area. The geoid 
model is given by 


М = МХ,Ү,2) +АМ +95 + п 


[5 


where 
S : the signal 
n :the noise 
N,(X,Y,Z) + AN : the system function from the surface models at control 


points where both h and H are known 
Then, for a given N,(X,Y,Z) the method of collocation can be used to determine Sq 
at these control points and to predict 5, at other points in the local area. At any point 


in the local area the value of h can be determined by using differential GPS 
measurements between a known point and any other point by the equation 


h = h + Ah 
The value of H at any point is given by 
H=h - N,(X,Y,Z)-AN-S-n 
The value of H at any point in the local area, which depends on the value of h at 
the control points and the differential GPS measurements, can be predicted with an 
accuracy of tn. The general form of the observation equation in the method of 


collocation is [Jeyapalan, 1977]: 


x = AeX+S +n +08 
а q р 


where 
x: the vector of the observation (x = Ah - N¿(X,Y,Z) - H) 
А :a given rectangular coefficient 
X : the vector of the systematic parameters (họ, a, b, c, d, or họ, a’, b, c, 
dee) 
S, :a signal vector at q observation points 
n, : a vector of measuring errors, noise at q points 
5, : a signal vector at p unknown stations 
e :indicates matrix multiplication 
: the null matrix 
If 


Z =S +n 
q q q 
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Then 
х= АеХ + 7, + 0е5 
9 9 р 
then 


peu eA AG. x 


and 


5 


Cae eee (xu А.Х) 
P Pq q 


where the variance-covariance matrix C of the 2, апа 5, vectors is 


oe Se 
ра 4 
The essence of this method is that by some means a covariance matrix can 


be assigned to the signal. For noise it will be possible to assign a diagonal weight 


matrix. 
5, are the values of the signal at the interpolated stations. Suppose there are 


q observations (and values of 99^. P interpolated values of 5, and m model 


parameters. The covariance matrices are the following [Bomford, 1980]: 


(1) Co the expected covariance between the observed x's for all pairs of 


the q observations. Itis aq x q matrix. 
(ii) Cio the expected covariance between the signals for all pairs of the q 


observations. It is aq x q matrix. 
(111) Cng, the expected value of n at each station. Itis a q x q diagonal 


matrix. 
(iv) Cog the expected covariance between all pairs of mixed observed 


and interpolated signals. It is a p x q matrix. 
(у) С,» the expected covariance between the signals at pairs of 


interpolated stations. It is a p x p matrix. 


The variance of the noise, Cua can be estimated in the usual way according to 


the circumstances at each station, different types of instrument, etc. The noise at 
different points is independent; hence 
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where 6, is the standard error of the noise. 
The expected covariance between the observation x's, C, can be obtained 
by 


С=С +С 
4 па Sq 


because the signal is small. 
The C,, can be be computed from a simple function whose paramcters can 


be determined by using empirical data. The covariance function is positive and 
definite. These two characteristics are found in many functions. Functions 
commonly used for covariance may be constant, sinusoidal, Gaussian, 
exponential, exponential cosine, and exponential sine and cosine. In this thesis the 
sinusoidal function was used. The function is given by 


C e = B sin(kr) 
CI) C a + B sin(kr) 


where с is the standard deviation of the control stations and B, k are 
coefficients to be determined. 

The parameters can be determined by least squares, and residuals at each 
point can be computed. From the residual, the covariance between points at a 
distance (r) can be computed by 


2. SAM 
Cn == =] 


where V, 15 the residual at the center, V, is the residual at a point which is at a 


distance r from the center and n is the number of points. There are several 
methods of determining C(r) of which the concentric circle approach is the 
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simplest. The coefficients B and k can then be computed from the computed 
covariance. 

The expected covariance between the signals at pairs of interpolated stations 
E can be obtained by substituting the distance between the interpolated stations 


and control stations into the covariance function. 


The expected covariance between all points of mixed observed and 
interpolated signals C,, can be obtained from the covariance function for each 


distance from the interpolated stations to the control stations. 
In this thesis two cases were studied in solving for the geoid height model. 


Case 1. C, = 1 and Era = 0 


Case 2. C, + 1 and Cx z 0 


I start with assuming 


С = Рр! = 
9 


С = 0 
Pq 


where P is the weight matrix and I is the unit matrix. 
Then 
X -(A PA)! A'Px 


and 
S, =0eP (x- AX,)=0 


Using X to compute residuals, v, and then estimating Ca C e and P, it is then 
found that 
OIE APTA EX 
1 4 4 


апа 
= и 
5, = Cua C (x-A X) 
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IV. DATA COLLECTION AND INSTRUMENTATION 


A. PLANNING 
Nine temporary bench marks were established on the NPS campus (Figure 4). 
Bench marks GH7 and GH8, designated as check marks, were established in the 


center of the levelling loop. H of bench marks was obtained by differential 
levelling. Ah was measured with a pair of GPS receivers. 
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Figure 4. The temporary bench marks on the NPS campus. 
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Eight permanent bench marks were recovered in the study area (Figure 5) 
[Vertical Control Data, 1961]. Bench mark GWM 27, designated as a check mark, 
was roughly in the center of the survey area for the permanent bench marks. H was 
obtained from the published elevations. Ah was also measured with a pair of GPS 
receivers. 
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Figure 5. The permanent bench marks in the Monterey Bay area. 


B. DATA COLLECTION 
1. Obtaining H 
a. Precise Levelling 

Precise levelling was run on the NPS campus on February 19, 1988, 
and March 24, 1988. The elevation of station TREE was assumed to be zero 
above the MSL, so the orthometric height differences, AH, for each temporary 
bench mark could be observed. The differences of the relative elevations, AH, 
for the forward and the backward sights on the NPS campus are listed in Table 1. 


Table 1. AH FOR THE NPS CAMPUS 


AH (m) 





To avoid an obstruction near bench mark D 697, which is 10 ft north 
of a 12-inch cedar tree, the temporary bench mark TEMP 1 was established 
nearby. Offset levelling was run on April 8, 1988. The elevation of TEMP 1 
was offset by -0.401 m from the elevation of D 697. 

F 813 which was set in the west end of the south abutment of a steel 
bridge over the Salinas River was offset to temporary bench mark TEMP 2. The 
offset levelling was run on April 9, 1988. The elevation of TEMP 2 was offset 
by -2.218 m from the elevation of F 813. 
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b. Published Elevation 

H values for the larger of our two study areas in and near Monterey, 
California, areas were obtained from the published elevations printed by U.S. 
Department of Commerce Coast and Geodetic Survey Washington D.C. [1961]. 
The geodetic datum used in this publication was the NGVD 29. First-order spirit 
levelling has extended this datum over most of the continent. Although first- 
order lines may be 300 km apart in some western areas, most points in the 
country are no more than 50 km from an estimated first-order bench mark. A 
readjustment of the this network is the NAVD 88. This new adjustment will be 
made to the geopotential surface rather than to sea level, and it will place the 
existing vertical data in a form that makes it most consistent and accessible to the 
user. Changes to older published elevations are not expected to exceed 15 
decimeters [NASA, 1978]. Table 2 lists the orthometric heights of the permanent 
bench marks we used. 


Table 2. H FOR PERMANENT BENCH MARKS UESD IN THIS 
STUDY 





2. Obtaining Ah 
a. Satellite Observation Plan 
Due to the positions of the satellite orbits during this study, the 
observing window was between 2100 and 0200 hours Pacific Standard Time 
(PST=UTC+8 hours.). The time period was between February 2, 1988, and May 
12, 1988. 
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b. Satellite selection 

The same five satellites (SV) (6, 9, 11, 12 and 13) were used for all 
observations. 

c. Position Dilution of Precision (PDOP) 

The accuracy to which positions are determined using GPS depends 
on two factors: (i) satellite configuration geometry, and (ii) measurement 
accuracy. GPS measurement accuracy represents the combined effect of 
ephemeris uncertainties, propagation errors, clock and timing errors, and 
receiver noise. 

The effect of satellite configuration geometry is expressed by the 
dilution of precision (DOP) factor, which is the ratio of the positioning accuracy 
to the measurement accuracy [Wells, 1987]. 


6 = DOP eo) 


where 
Со : 15 the measurement accuracy (standard deviation) 
© :1s the positioning accuracy (standard deviation in one coordinate) 
The value of GDOP itself is a composite measure that reflects the 
influence of satellite geometry on the combined accuracy of the estimation of 
observation time (user clock offset) and receiver position [Milliken, 1986]. 


GDOP/- PDOP + TDOP 


TDOP is the Time Dilution Of Precision, the error in the clock of the receiver 
bias multiplied by the velocity of light. The four best satellites selected by the 
receiver are those with the lowest GDOP. Trimble recommends that the rapidly 
changing PDOP provides better geometry for phase differencing techniques. 
Low or constant PDOP provides weaker solutions. The 4000SX receiver does 
not record GDOP, but it does record PDOP every five minutes. PDOP was 
about 4.5 for all GPS observations. 
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C. INSTRUMENTATION 
1. Levelling 

A Zeiss model Ni-2 level (Ser. # 82377) was used at the temporary 
bench marks for the third-order, class I levelling. It has a 32-diameter 
magnification, produces an erect image, and has stadia constants of 333 or 100. 
A Peg test was performed before levelling. The level error, or c-value was also 
checked before the beginning the levelling. The c-value was -0.006 mm/m which 
was less than +0.05 mm/m, so it was not necessary to adjust the level [Bodnar, 
1975]. The level contains a bubble tube to permit positioning parallel to the 
geoid. When properly set up at a point, the telescope is locked so that it will 
rotate through a 360° arc in a horizontal plane. With the level locked in position 
readings are made on two calibrated staffs held in upright positions ahead of and 
behind the instrument. The difference between readings is the difference in 
elevation between points. Dietzgen Model # 6450 metric rods were used in the 
levelling. These rods are graduated in centimeters. The actual reading is 
estimated to the nearest 1 mm. Rod levels were used to indicate when were 
vertical. 


2. GPS 
a. Trimble 4000SX Receiver 

A complete description of the Trimble 4000SX receiver is given by 
Trimble Navigation [Trimble, 1987a]. NPS operates three Trimble 4000SX GPS 
receivers. For this study I fixed one antenna to the roof of Building 224 on the 
NPS campus, and one was carried to the field. The 4000SX is capable of 
observing the C/A code, integrated Doppler, and carrier beat phases of up to five 
satellites simultaneously. Its ability to use the C/A code allows the receiver to be 
used as a Stand alone navigation system to determine positions using Doppler- 
smoothed pseudoranges and velocities [Ashjaee, 1985]. The receiver uses the 
C/A code in a time transfer mode to determine the offset and drift of its own 
clock and thus provide accurate time tags for the observations without the 
requirement of an external atomic clock or synchronization with the receiver at 
the other end of the baseline. The reference position (the geodetic coordinates of 
the antenna) and the practical options chosen must be entered into the receiver via 
the receiver key pad. The 4000SX requires 115 V AC power or 12 V DC 
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power supply. 115 V AC power was used for the NPS campus measurements. 
12 V DC, supplied by a car battery, was used in the field. 
b. Grid Personal Computer (PC) 

For precise relative positioning the 4000SX receiver transmits data 
through an RS-232 port to a microcomputer (Grid PC) for storage on floppy 
disks for post-processing. The Grid PC uses either 115 V AC or 12 V DC. On 
the NPS campus a 115 V AC power supply was used, so measurements could be 
made for four hours. The 12-volt battery in the Grid PC lasts about 120 
minutes, SO measurements were taken for only 100 minutes in the field. 

c. Antennas 

Multipath-resistant Trimble microstrip antennas were installed on 
Building 224 on the NPS campus and over the various bench marks. The antenna 
heights from the center of bench marks to the edge of the antenna's ground plane 
were measured before and after GPS observations. The field antenna was 
mounted on a tripod with a tribach and optical plummet for centering and 
levelling the antenna. Arrows on the antennas were oriented to the north at both 
stations using a magnetic compass. 

d. Meteorological Instruments 

A barometer and a sling psychrometer were used to measure 

atmospheric pressure, relative humidity and air temperature at each field site. 
e. The steel tape 

A three-meter steel tape was used to measure the antenna height. 

This tape can be read to 1 cm. Readings are estimated to 1 mm. 


D. SOFTWARE 

A complete description of the Trimble-supplied Trimvec software is to be 
found in Trimble Navigation [Trimble, 1987b]. The software provides data 
logging, baseline computation and datum transformation programs. These 
operate with an IBM compatible personal computers. A printer and a hard disk 
are used with the planning and processing programs. The following programs 
are used on 3-1/2 inch micro-floppies: 

1. Data Logger on Disk 1 

In relative positioning, the receiver is controlled from the Grid PC by 

version D of Trimble's Gridlog5.bat program. Each observation session was 
initialized to log data when a minimum of four satellites were 15° above the 
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antenna's horizon. Five satellites were designated for each observing session. 
The observables and receiver clock parameters were logged to a floppy disk 
every 15 seconds and the C/A code-determined antenna position and PDOP every 
five minutes. The GPS navigation message was logged to a separate file at the 
beginning of the session. 

2. Post-processor on Disk 2 

The data was processed using the Trimvec Trim640 program, Revision 
AB. Trim640 is a relative positioning, post-processing program that provides 
triple and double difference solutions for two sets of the 4000SX carrier phase 
data logged simultaneously at two stations. Trim640 adopts the best C/A code 
position during the data loading. Only the broadcast ephemeris can be used to 
compute fixed orbit satellite positions. Trim640 limits processing to 700 
epoches, so the first 700 epoches for each observing session are used. The 
antenna on the roof of Building 224 on the NPS campus was used as reference 
station and its coordinates were kept fixed. Permanent bench marks were chosen 
as Trim640 reference stations for the off campus area because the observing 
sessions at them were shorter than the observing session at Building 224 on the 
NPS campus. 

3. Geoid.exe on Disk 3 

This program computes global geoid height, Ny(X,Y,Z), at any point 
with an accuracy of a few meters [Trimble, 1987b]. The global geoid heights, 
Ny(X,Y,Z), are based on Rapp's 1978 360 x 360 model, an harmonic expansion 
of the geopotential referred to GRS 80. Heights are suitable for showing relative 
shape and trends in geoid. Global geoid height for the Monterey Bay area from 
36° 34° N to 36° 49' N latitude, and from 121? 34' W to 121° 55' W longitude is 
given in Figure 6. A Fortran program GPSCON (Appendix A) was used for the 
contoured plot. 

A second program also named Geoid.exe was obtained from the NGS of 
NOAA's National Ocean Service Charting and Geodetic Services uses the same 
Rapps 1986 360 x 360 model. The global geoid height in the Monterey Bay area 
from the NGS program is shown in Figure 7. Because the NGS program rounds 
off to the nearest decimeters, there are some steps on the curves. The differences 
between Figure 6 and 7 may be attributed to differences between Rapp's 1978 
and 1986 models. 
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Figure 6. The global geoid height in the Monterey Bay area calculated 
using the Trimvec Geoid.exe program. 
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Figure 7. The global geoid height in the Monterey Bay area calculated 


using the NGS Geoid.exe program. 
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4. Satellite Visibility Program on Disk 3 


The program Stvis.exe is a system for generating visibility data for the 
GPS satellites. The program takes as input the observer's latitude, longitude, a 
mnemonic name for the location, data for calculations and optionally the time zone 
offset from Universal Time Coordinates (UTC). NPS in Monterey, California, was 
used as the reference station ( 36? 35' 56.1" N, 121? 52' 36.0" W, and altitude -19 
m). Satellite orbit data is contained in the file, Almanac.dat, which is produced by 
processing data collected from the actual satellites. The satellite visibility chart 
shows the tracks of the GPS satellites during the planned observation session 


(Figure 8). 





Figure 8. Sky plots of satellite tracks for the Monterey Bay area. 
Elevation angles are dotted concentric circles. Zenith is at the center. 
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V. DATA PROCESSIN 


A. ORTHOMETRIC HEIGHT ON NPS CAMPUS 

In calculating the most probable elevations for the temporary bench marks, 
station TREE was used as a reference, where the elevation was assumed to be O 
meters above MSL. The interlocking levelling circuit on the NPS campus is 
shown in Figure 9. 





Figure 9. The interlocking levelling circuit on the NPS campus. 


A constrained least squares adjustment was used to do the computation. A 
weight of 100 was assumed at the TREE and weights of 1 were used for all 
observations. The 13 observation equations are: 
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Pa Ж P ( 0.000 + V, ) 
P, (H, - Hga) = P, ( -0.808 + V, ) 
POLI Р ( 5.386 + V, ) 
P (H, - H,) =P,( 3.149 + У 
РОН = P, (11.535 + V 
P. ( H,- H, ) = Р, ( -1.016 + У, ) 
Р, (Н, - Н, ) - P,( 4214 * V, ) 
P ( H,- H,) =P,( 0489 + V, ) 

) 

) 


P, (Higep- HH, ) = P, ( -4.032 + V 
PCH H) =P 0.547 + У 
P ( H,- H,) =P ( 0.057 + V) 
РИН ан Б 0717 + у.) 
Se ORI Y) 


The observation equations in matrix form are [Wolf, 1980]: 
PAH= PL + PV 


The weight matrix P,4.,,, coefficient matrix A,3,,, observation matrix L,,,,, 
and residual matrix V,,,, are: 
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The program Lobs.basic [Jeyapalan, 1988] was used to find H for the temporary 
bench marks. The program (Appendix B) gives the H and V matrices: 
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0.00000 


0.00007 

0.000 0.00007 
-0.808 -0.00022 
4.579 -0.00044 
7.128 -0.00022 
He 59 769 У = | -0.00022 
8.246 0.00007 
4.032 -0.00019 
4.521 -0.00010 
8.445 -0.00010 
0.00022 

0.00023 


The standard deviation of the observation equations is 0.36 mm. 


B. GPS DATA PROCESSING 
1. Automatic Processing Mode 

Trimble recommends that an automatic processing mode be used if more 
than one hour of data from four or five GPS satellites is taken at each site for 
baselines lengths up to 30 km. 

The output files contain triple difference, double difference float and 
double difference fixed solutions. The double difference fixed solution provides 
the most accurate solution if the following conditions are met: 

a. Integer bias search indicates the ratio sum-of-squares must be greater than 
3.0; 

b. RMS (cycles) less than 0.05; and 

c. Difference between the float and fixed solutions is less than 10 cm, in any 
component of the baseline (X,Y,or Z). 

In the GPS data processing all the solutions from the bench marks 
satisfied these conditions except station J 697 for which RMS was 0.068. The 
distance from NPS Building 224 to J 697 is 33 km. Trimble recommends a 
baseline length greater than 30 km when the fixed solution does not meet the 
above conditions. The float solution should be used provided the RMS of fit is 
better than 0.08. Ah of double difference fixed solutions were used for the bench 
marks in this study. The double difference fixed solutions are shown in Table 3. 
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Table 3. DOUBLE DIFFERENCE FIXED SOLUTIONS 


Slope distance Coordinates difference between 
from Bldg. 224) RMS Ratio fixed and float solution 
cycles) 





Default parameters are assumed for automatic processing. These 
parameters have a minimum elevation mask of 15° for double differences. 
Station one coordinates were the best code positions during data logging. Actual 
surface meteorological values were used in the processing. 

For sets of data covering more than one hour, epoch increments of five 
or ten provide full precision. Automatic processing begins with several iterations 
of triple differences to establish a starting value for the double difference 
solution. A cycle slip fixer processes every epoch. 

Double difference solutions begin with several iterations of the float 
solution where biases, as well as baseline components, are computed. It is called 
the float solution, since the biases are allowed to float and are not constrained to 
be integers. Each iteration should indicate convergence toward zero for 
observed minus computed phases. The default elevation mask is 15° for 
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processing double differences. After completion of the float solution, 
correlations and biases are computed. The integer search algorithm begins to 
determine the proper integer values for the bias. After the biases are set to their 
integer values, the processing is repeated to compute the fixed solution. Ellipsoid 
height differences Ah for the temporary bench marks with respect to NPS 


Building 224, found using double difference fixed solutions is given in Table 4. 


Table 4. Ah AND WGS 84 COORDINATES FOR THE 
TEMPORARY BENCH MARKS. 


am | X - 


X (m) 
-2707298.015 
-2707385.942 
-2707502.139 
-2707528.391 


-2707367.536 
-2707202.581 
-2707252.929 
-2707406.088 
-2707329.452 


-4353450.286 
-4353407.659 
-4353548.202 
-4353750.407 
-4353813.483 
-4353799.531 


--4353587.524 


-4353553.901 
-4353759.387 


3781781 .432 
3781790.139 
3781549.491 
3781314.047 
3781306.605 
3781493.061 
3781666.956 
3781589.428 
3781428.906 





Ah for the permanent bench marks with respect to NPS Building 224 
using double difference fixed solutions are given in Table 5. 
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Table 5. Ah AND WGS 84 COORDINATES FOR THE 
PERMANENT BENCH MARKS 


















































2.2398 -2696890.306 | -4350931.132 | 3792065.393 
-8.4275 -2708436.436 | -4351973.362| 3782674.213 

1.6335 -2692137.111 | -4360870.528 | 3784094.277 
-1.3569 -2689245.909 | -4360679.825 | 3786312.651 
15.1728 -2685153.665 | -4357222.820 | 379317541171 
71.3023 -2679276.491| -4356508.580| 3798216.777 
-6.7070 -2695613.779 | -4350454.843 | 3793463.222 


49.1051 -2692408.721 | -4360427.654 | 3784433.846 


2. Obtaining Geoid Height 


Geoid heights of bench marks were obtained by using the Geoid.exe 
program from NGS or Trimvec (Table 6). 


Table 6. GEOID HEIGHTS 
-34.30428 
-34.30671 
-34.31508 
-34.32072 
-34.31605 
-34.30745 
-34.30532 
-34.31140 
-34.31247 
-33.83303 
-34.31918 
-33.86473 
-33.76630 
-33.58265 
-33.42047 
-33.77899 
-33.86547 
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3. Determination of Local Geoid Model 
The method of collocation was used to solve the local geoid model. The 
procedure is described as follows : 
a. Reference Station Selection 
Station TREE was selected as a reference station on the NPS campus. 
Coordinate differences AX, AY and AZ in WGS 84, for the temporary bench 
mark with respect to the TREE were computed. These are shown in Table 7. 


Table 7. AX, AY AND AZ ON THE NPS CAMPUS 


AX (m) 
0.000 
-87.927 
-204.124 
-230.376 
-69.521 
95.434 
45.085 
-108.073 
-31.437 


0.000 
42.627 
-97.916 
-300.121 
-363.197 
-349.245 
-137.238 
-103.615 
-309.101 


0.000 

8.707 
-231.941 
-467.385 
-474,827 
-288.371 
-114.476 
-192.004 
-352.526 





Station K 152 was selected as a reference station in the large off 
campus area. The coordinate differences, AX, AY and AZ between the 
permanent bench marks and K 152 were computed in WGS 84 (Table 8). 
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Table 8. AX, AY AND AZ IN THE OFF CAMPUS AREA 


AX (m) 


0.000 
-11546.130 
4753.195 


7644.397 
11736.641 
17613.815 

1276.527 

4481.585 





0.000 
-1042.230 
-9939.396 
-9748.693 
-6291.688 
-5577.448 

476.289 
-9496.522 


0.000 
-9391.180 
-7971.116 
-5752.742 

1109.778 
6151.384 
1397.829 
-7631.547 


b. Determine the Local Geoid Model 
I start with assuming the signal, S, equals zero and the weight matrix 
P, equals the unit matrix. Seven control marks in both areas led to seven 
observation equations. The coordinate differences between the reference marks, 
which are TREE on the NPS campus and K 152 in the off campus area, are the 
coefficients of the parameters. The observation equation is given by 


AN = H +N (X,Y,Z) - Ah 


where 


H : the orthometric height, from levelling 
N,(X,Y,Z) : the global geoid height, obtained form NGS or Trimvec 


Ah : the ellipsoid height difference with respect to the reference station 
The local variation of the geoid height, AN, is also given by 


AN = hy + combination of coordinate differences 


hs, the ellipsoid height, includes a constant correction to N,(X,Y,Z) at the 
reference station. AN for temporary bench marks, N,(X,Y,Z) from the NGS 
(designated as NAN) and N,(X,Y,Z) from the Trimvec (designated as TAN) are 


given in Table 9. 
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Table 9. AN ON THE NPS CAMPUS 


TAN (m) 


-26.3468 -25.9511 
-26.3266 -25.9333 


-26.3444 -25.9595 


-26.3338 -25.9545 
-26.3435 -25.9595 
-26.3441 -25.9515 
-26.3283 -25.9336 





AN for permanent bench marks with N,(X,Y,Z) from the NGS (designated as 
NAN) and N,(X,Y,Z) from the Trimvec (designated as TAN) are given in Table 


10. 


Table 10. AN IN THE OFF CAMPUS AREA 
-19.3748 -19.1078 
-20.4855 -20.1047 
-19.5555 -19.5202 


-19.3091 -19.2754 
-18.7158 -18.6984 
-19.5413 -19.6618 
-19.3100 -19.0890 





A least squares adjustment program Lobs.basic was used to determine 
the geoid model which has the smallest standard deviation. This was done by 
using the different combinations of the coordinate differences. GPS data from 
the off campus area were studied to determine the local geoid model. N,(X,Y,Z) 


from the Trimvec program were used to compute AN. 
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Four-parameter combinations and their standard deviations for the off 
campus area are shown in Table 11. 
Table 11. FOUR PARAMETERS 


Standard deviation (m 
0.45 


0.52 
0.34 
0.47 
0.52 
0.34 
0.47 
0.31 
0.37 
0.52 
0.34 
0.46 
0.43 
0.58 
0.28 
0.37 
0.32 
0.49 
h, + aAY* + bAX? + cAXAY 0.20 * 


* The smallest standard deviation. 


Parameters 
hy + aAX + bAY + cAZ 


h, + aAX + bAY  cAY? 

hy + aAX + bAY + cAX? 

hy + aAX + DAY + cAXAY 
h, + aAX + bAZ + cAY” 

h, + aAX + DAZ + cAX? 

ho + aAX + bAZ + cAXAY 
h, + aAX + bAY? + cAX° 
D. + aAX  bAX?  cAXAY 
h, + aAY + bAZ + cAY° 

hy + aAY + bAZ + cAX? 

hy + aAY + bAZ + cAXAY 
h, + aAY + BAY’ + cAX” 
h, + aAY+ bAY’ + cAXAY 
h, + aAY+ bAX’ + cAXAY 
hy + aAZ + bAY? + cAX? 
h, + aAZ + bAX” + cAXAY 


h, + aAZ + bAY? + cAXAY 
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Five-parameter combinations and their standard deviations for the off campus 
area are given in Table 12. 


Table 12. FIVE PARAMETERS 


Parameters Standard deviation (m) 
h, T HAX + bay + CAZ + daR’ 


hy Pa F DAY F cAZ + da 
iy к аАХ + БАҮ + AZ + dAXAY 


hy FAN + BAZ e ca дАХДҮ 
hy TRAY + DAZ + cAX? + dAXAY 
hy aA + ЪА7 + сАХТ + dAXAY 



















* The smallest standard deviation. 


Six-parameter combinations and their standard deviations for the off campus area 
are given in Table 13. 


Table 13. SIX PARAMETERS 


Parameters 
h, + aAX + BAY + cAZ + dAX’ + eA Y 
h, + aAX + BAY + cAZ + dAX* + eAXAY 
h, + aAX + bAY + cAX? + dAY? + eAXAY 
h, + aAX + bAZ +cAX” + dAY” + eAXAY 
h, + aAY + bAZ + cAX? + dAY? + eAXAY 


* The smallest standard deviation. 
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For seven observation equations and seven parameters the degree of freedom 
equals zero, so the standard deviation is undefined. Seven-parameter 
combinations and their standard deviations for the off campus area are given in 
Table 14. 


Table 14. SEVEN PARAMETERS 


Standard deviation (m 
h, + aAX + bAY + cAZ + dAX’ + елү? + ҒАХДҮ 


The best standard error of the four-parameter geoid model is greater 
than for five-parameter or six-parameter geoid models. The seven-parameter 
geoid model did not converge, so five-parameter and six-parameter geoid models 
were chosen for this study. The five-parameter geoid model is given by 


N = N,(X,Y,Z) + h, + aAY +bAX” + CAY” + dAXAY 


and the six-parameter geoid height is given by 


N - NQOLY,Z) * h, - a'AX ^ b'AZ + с'АХ°* + d'AY” + e'AXAY 


c. Evaluating the Geoid Model Using Check Marks 
In order to evaluate the accuracy of the geoid model, the geoid height 
model can be rearranged to calculate the orthometric height of check marks 
which are GH7 and GH8 on the NPS campus and GWM 27 in the off campus 
area. 


Н = hy+ Ah - N,(X,Y,Z) + aAY + bAX? + cAY? + dAXAY 


OT 
Н = hy+ Ah - N,(X,Y,Z) + a'AX + b'AZ + c'AX? + d'AY? + e'AXAY 


The ho, a, b, c and d of the five-parameter geoid model, and hy, a', b', c', d' and 


e' of the six-parameter geoid model, can be solved by least squares adjustments. 
The results for hp, a, b, c and d for the five-parameter model using the seven 
control marks in the off campus area are given in Table 15. The results for ho, 
а, б, с and d are shown in the NGS column, where the N,(X,Y,Z) was from the 
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NGS Geoid.exe program. The results for h,, a, b, c, d and standard deviation, o, 
are shown in the Trimvec column, where the N,(X,Y,Z) was from the Trimvec 


Geoid.exe program. 


Table 15. hj, a, b, c, d; 0 FOR SEVEN CONTROL MARKS 


-19.26895 -19.00915 


-2.344959E-04 -2.736598Е-04 


-8.523017Е-09 -8.528133E-09 
-3.574314E-08 -3.905052E-08 
-2.211550E-08 -1.728313E-08 
0.130518 0.095723 





The results of ho, a, b, c, d and o of the five-parameter geoid model 


using the six control marks (excluding B 21) in the off campus area are given in 
Table 16. 


Table 16. hy, a, b,c, d,o FOR SIX CONTROL MARKS 


NGS (m) 


-19.25415 -19.01411 


-2.752543E-04 


-7.735253Е-09 

-3.777950Е-08 

-1.798617Е-08 
0.171919 


-2.600850Е-04 
-8.792540Е-09 
-3.836067E-08 
-1.867193E-08 
0.133516 





The results of h,, a, b, c, d and o of the five-parameter geoid model 


using the six control marks (excluding J 697) in the off campus area are given in 
Table 17. 
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Table 17. h a, b, c,d, © FOR SIX CONTROL MARKS 


NGS (m 


-19.265355 -19.03487 


-2.503395E-04 -1.977533E-04 


-8.731945E-09 -7.526410E-09 

-3.711466E-08 -3.246532E-08 

-2.179058E-08 -1.882791E-08 
0.1841317 0.1206984 





The results of h,, a, b, c and d of the five-parameter geoid model using 
the five control marks (excluding B 21 and J 697) in the off campus area are 
given in Table 18. The standard deviation is undefined since degree of freedom 
equals zero. 


Table 18. h, a, b, c, d FOR FIVE CONTROL MARKS 
-19.37497 -19.10791 
1.032352E-04 3.397465E-05 
8.119969E-09 3.521564E-09 
7.377821Е-09 -3.339665Е-09 
1.338776Е-09 -3.667083E-09 





The results of họ, a, b, c, d and 6 of the five-parameter geoid model 


using the seven control marks on the NPS campus are given in Table 19. 


Table 19. hj, a, b, c, d; o FOR SEVEN CONTROL MARKS 


-26.33543 -25.94105 
2.976507E-06 1.192093E-07 


-5.419133E-08 -1.763692E-07 
-3.601599E-08 -8.489588E-08 
6.309711E-08 -3.702007E-08 
0.013773 0.014612 





The results of h,, a', b', c', d', e and o of the six-parameter geoid model 


using the seven control marks in the off campus area are shown in Table 20. 


Table 20. h,, a', b', c', d', e' c FOR SEVEN CONTROL MARKS 


NGS (m) 


-19.27388 -19.01660 


2.083 182E-04 1.8423 80E-04 
-2.561510Е-04 -2.520234Е-04 


-7.731387E-09 -8.335974E-09 

-3.767036E-08 -3.960304E-08 

-1.249646E-08 -1.521767Е-08 
0.137149 0.123960 





The results of hp, a, b', c', d' and e' of the six-parameter geoid model 


using the six control marks (excluding B 21) in the off campus area are given in 
Table 21. The standard deviation is undefined since degree of freedom equals 
Zero. 


Table 21. h,, a', b', c', d', e FOR SIX CONTROL MARKS 


NGS (m 


-19.37489 -19.10785 
2.5764 1 Е-04 2.288222Е-04 


-1.691282Е-04 -1.734756Е-04 
-1.083072Е-08 -1.113654Е-08 
-2.817797Е-08 -3.102741Е-08 
-5.824404E-09 -9.185896E-09 





The results of ho, a’, b‘, c', d' and e' of the six-parameter geoid model 


using the six control marks (excluding J 697) in the off campus area are given in 
Table 22. The standard deviation is undefined since degree of freedom equals 
zero. 
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Table 22. hy, a’, b', c', d', e' FOR SIX CONTROL MARKS 


NGS (m) 


-19.37466 -19.10774 


1.065433E-04 9.238720E-05 


-4.604459E-05 -6.231666E-05 
-2.012712E-09 -3.170499E-09 
-1.141234Е-08 -1.586523Е-08 
-2.528395E-09 -6.206392E-09 





The results of ho, a, b’, c', d', e апас of the six-parameter geoid model 


using the seven control marks on the NPS campus are given in Table 23. 


Table 23. hy, a’, b’, c’, d’, e’, o FOR SEVEN CONTROL MARKS 


-26.33389 -25.93820 
1.645088E-05 4.556775E-05 
4.225970Е-05 6.169081Е-05 


6.856863Е-08 6.391201Е-08 
6.007031Е-08 5.878974Е-08 
1.781154Е-07 1.767185Е-07 
0.018859 0.018872 





Fortran program GPSDIS (Appendix C) was used to calculate the 
orthometric heights of the check marks. The results are discussed in the Chapter 
VII. 

d. Obtaining the Weight Matrix 

Since the standard error of the geoid model is about + 2 cm, which is 
the same or better than GPS Ah accuracy, it is not actually necessary to do 
further computation. Nevertheless, the weight matrices, P, were computed. 

The five-parameter geoid model used for the NPS campus and the 
six-parameter geoid model used for the off campus area were employed in 
obtaining P. The covariance function, C(r), of the control marks (p. 18) can be 
obtained by using the least squares residuals. The coefficients B and k of C(r) 
can be obtained by solving simultaneously the equations for C(r) with different 
distances from the centers. 
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(1) P. for NPS Campus. The residuals of the five-parameter geoid 
model of the control marks on the NPS campus are given in Table 24. The 
global geoid heights were from the NGS. 


Table 24. RESIDUALS FOR FIVE-PARAMETER MODEL 


# of residual Residual (m) 


1.137352E-02 
-9,420395E-03 
7.339478E-03 


-4.278183E-03 
3.572464E-03 
6.446839E-04 

-8.714676E-03 





For r= 593 m 


рату РУ FV OVE + VV 
Ст) E E. A 2—2 = -2.194272E-03 


For Г, = 319 т 


МУ Уа баға дары е АЙ E274 VA V АА ТАУ, ay 


C(r,)= 3 2 — 3.945813E-6 


Thus the equations of covariance can be written as 
С(т,) = 0.0138+В sin(kr, ) 


and 
C(r,) = 0.0138 + B sin(kr,) 


Approximate solutions for coefficients B and k can be found by using the 
Maclaurin series expansion. The solved covariance function is given by 


C(r) = 0.0138 - 0.0160 sin(1.7670 г) 
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C. can be obtained by using the C(r) with the distances between the control 
marks. Fortran program DISTCO (Appendix D) was used to compute the 
distances between the control marks, as well as с? and, hence, the weight matrix, 
p Cay The Matlab program on the NPS mainframe computer was used to find 
the inverse of C. (P= CEDI P for the five-parameter geoid model for the NPS 
campus is given by 


162 ZEROS 


ZEROS 100 


(2) P for off Campus Area. The residuals of the six-parameter geoid 
model of the control marks in the off campus area are given in Table 25. The 
global geoid heights were calculated using the Trimvec program. 


Table 25. RESIDUALS FOR SIX-PARAMETER MODEL 


# of residual Residual (m) 


9.120178E-02 
-0.773254E-03 
6.391526E-03 


1.983643E-04 
-2.779961E-02 
1.685524E-02 
-7.651711E-02 





For г = 13994 т 


УУ + УУ + УУ + УУ. +УУ + УУ 
Cua) E a = -6 659883E-04 
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For г, = 19230 т 


У ММУ. 


C(r,) = 5 


= 7.912463E-04 


Thus the equations of covariance can be written as 


C(r,) = 0.0124 + В sin(kr ) 


C(r,) = 0.0124 + В sin(kr,) 


Approximate solutions for coefficients B and k can be found by using a 
Maclaurin series expansion. The solved covariance function is given by 


C(r) = 0.0124 - 0.1329 sin(2.8278 r) 


Similarly, C, can be obtained by using the C(r) with the distances between the 


control marks. Fortran program DISTCO (Appendix D) was used to compute 
the distances between the control marks and C,, to obtain the weight matrix (P = 
E The Matlab program on the NPS mainframe computer was used to find the 
inverse of C, (P 2 C,!). P for the six-parameter geoid model of the off campus 
area is given by 


25 ZEROS 
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VI. EVALUATION OF GPS RECEIVER ERRORS 


Factors affecting the GPS measured ellipsoid height differences, Ah, 
include: (1) comparing GPS differencing solutions, (2) standard error of GPS 
observations, (3) corrections for surface meteorological values, and (4) 
observation durations for GPS. 


A. COMPARING GPS DIFFERENCING SOLUTIONS 
1. Ellipsoid Height Differences Tested at a Fixed Position 
In order to determine the best solution for Ah using GPS, the Ah data 
from a Trimble 4000SX GPS receiver were collected in the K parking lot on the 
NPS campus from 0633 to 1007 UTC on February 19, 1988. The design of the 
experiments is shown in the Figure 10. 


12-foot ladder 


Adjustable 
height 


Trimble 40005Х 
C 


Plumb bob 


Fixed point on the ground 





Figure 10. Illustration of the Ah test of Trimble 4000SX GPS 
receiver. 


50 


A plumb bob was used to point to a fixed position (36° 35' 56" N, 121° 
52' 36" W) on the ground. The observations were taken by changing the antenna 
height over a distance of about 60 cm, which was measured with a tape every 
half hour. The ellipsoid height differences, Ah, and the orthometric height 
differences, AH, should agree with the change of antenna height for a fixed 
position. Triple difference solutions, Ahypy, are given in Table 26. Table also 
gives the differences between the Ah, at succeeding observation times, DAh..gy; 
double difference float solutions, Ahp, y; the differences between the Ahp; y at 
succeeding observation times, DAh,, 7; and double difference fixed solutions, 


Аһу, апа the differences between the Ah,,, at succeeding observation times, 
DAhgyx 


Table 26. Ah ne an UR a AURA HEIGHTS 


E EL 
(19-2-1988) 





The differences between the GPS DAh and the differences of the 
measured ellipsoid height differences, DAh,,,,, are given in Table 27. 
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Table 27. THE DIFFERENCE OF Ah 


ANTENNA gr A |DAhyg, - DAhpp, | DANY, - БАһат|РАһ,у рл - РАр 
HEIGHT (m) m (m m) 





The results show that the double difference fixed solution is the best. 
The difference between the measured ellipsoid height difference and the GPS 
ellipsoid height difference was about 1 mm. The last two results were not good 
because they had multipath effects and imaging effects from the ground (antenna 
heights 1.074 m and 0.484 m), so the solutions of AH compared to the true 
differences were large. This indicates that the antenna should be at least 1 m 
above the ground. 

2. Ellipsoid Height Differences Tested at Different Positions 

Since the geoid on the NPS campus has a flat geoid slope, the differences 
between Ah, DAh, and the orthometric height difference, AH, of the temporary 
bench marks can be neglected. The Ah and DAh from the GPS observations are 
given in Table 28. The differences between the GPS DAh and the AH are given 
in Table 29. 
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Table 28. Ah AND DAh ON THE NPS CAMPUS 


mark m) (m) (m) (m (m) (m) 





Table 29. DIFFERENCES BETWEEN DAh AND AH 


mark (m) (m) (m) (m) (m) 
0.000 
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It is clear that the double difference fixed solutions give the best solution 
in a flat geoid slope area. 


B. STANDARD ERROR OF GPS OBSERVATIONS 

To determine the accuracy of GPS, measurements were taken at bench mark 
TREE on the NPS campus. Three measurements, each of about three hours 
duration, were taken for this analysis. Double difference fixed solutions were 
used to analyze the data. The double difference fixed solutions for these 
measurements are given in Table 30. 


Table 30. Ah... AT STATION TREE 


Feb. 5, 88 0730 - 1025 


Feb. 7, 88 0721 - 1004 
Feb. 11, 88 0705 - 0959 





The standard error of the mean of Ah (o / Vn) is about + 2 cm at TREE and 
the standard error (0) of the single measurement is about + 4 cm for three hours 
of observation. 

Results obtained by Strange in 1985 for a project southeast of Phoenix, 
Arizona, show that using multiple GPS observations at the same point lead to 
uncertainties typically less than 2 cm in Z-direction over 20-km distances 
[Zilkoski, 1988]. 

Thus the accuracy of the ellipsoid differences obtainable by GPS 
measurements is expected to be about + 2 cm for about three hours of 
observation. 


C. CORRECTIONS FOR SURFACE METEOROLOGICAL 
VALUES 


The meteorological correction is a function of the atmospheric refractivity 
as computed by the surface meteorological values of pressure, temperature, 
humidity and the elevation angle of the satellite. Larger corrections are required 
for low elevation angles, as the signal travels a longer path through the 
troposphere. A modified Hopfield troposphere model [Fell, 1975] is used for 
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automatic processing of Trim640 software. The model corrects for at least 90% 
of the tropospheric delay [Remondi, 1984]. Default surface meteorological 
values for automatic processing are 1010 millibars, 20° C and 50 % relative 
humidity. These parameters can be reset before the processing. Trim640 allows 
only one pressure, temperature and humidity entry for each site per session. 
Ellipsoid height differences, Ah,,,, were calculated for four marks on the NPS 
campus and two off campus using the default meteorological parameters and 
using the true values. There are given in Table 31 along with differences found 
between the two methods. 


Table 31. COMPARISON OF SURFACE METEOROLOGICAL 
CORRECTION 


Default Ah, (m) | True Аһ. (т) Difference (mm) 


-3.7777 -3.7766 


-0.6376 -0.6383 


0.9070 0.9055 
-3.8524 -3.8509 
2.2398 2.2430 
-71.3179 -71.3023 





The difference is about + 2 mm on the NPS campus and about + 2 cm in the 
Monterey Bay area. The Trim640 program contains no ionospheric model. For 
first-order relative positioning (1 part in 10°) and a baseline shorter than 50 km 
ionospheric delay can be neglected [Remondi, 1984]. In this study surface 
meteorological corrections were made during data processing to get the best Ah 
solutions. 


D. OBSERVATION DURATIONS FOR GPS 

To find out how long observations must be made to obtain the best solutions 
in the field when batteries are used for power, the GPS data from S 812, J 697 
and K 152 were segmented into averaging periods of length nAt, where At = 10 
minutes and n = 1, 2, 3, ...,10. Default meteorological values were used for this 
processing. Ah,,, for various observation durations for the stations is given in 


Table 32 and plotted in Figure 11 to Figure 13. 


DD 


Table 32. y AS A FUNCTION OF OBSERVATION DURATION 


S 812 (m) J 697 (m) K 152 (m) 


7529] 
70.9697 
71.2185 
71.4514 
71.4762 
71.5343 
71.5309 
71.3407 
71.3179 


Ah (m) 1.60 


40 60 
Time (min) 





Figure 11. Ah vs. observation durations at station S 812. 
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40 60 
Time (min) 


40 60 
Time (min) 





Figure 13. Ah vs. observation durations at station K 152. 


It is clear that Ah was affected by the observation duration. Ah varied quite 

a lot in the first forty minutes. After about forty minutes, Ah tended to close to 

the mean value of the observations. After sixty minutes, Ah didn't vary much. 

Trimble recommends, to ensure sufficient quality data and to obtain first-order 

results on single observations, at least one hour should be taken for a four or five 
57 


satellite session. For baseline lengths up to 30 km the conditions for the fixed 
double solution of automatic processing will be met if data are collected for 
ninety minutes on four or five satellites [Trimble, 1987b]. The results shows that 
the best solutions can be improved with longer observation times. Normally, at 
least sixty minutes of observations are required for vertical control. 
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VII. RESULTS AND DISCUSSION 


A. EVALUATATION OF GEOID MODEL 


To find the accuracy of the geoid model the orthometric heights of the 
check marks from GPS, Hops, were calculated, and these were compared to the 
orthometric height of the check marks from the levelling, H; pg, ; m. 

I started with the assumptions that C, = 1 and C,, = 0. Five-parameter and 


six-parameter geoid models were used to calculate the orthometric height of the 
check marks. GWM 27 is a check mark in the off campus area. GH7 and GH8 
are the check marks on the NPS campus. N,(X,Y,Z) was calculated using both 
the NGS and the Trimvec programs. 
1. Hops from Five-Parameter Geoid Model 
a. Hops in the off Campus Area 

Hgps in GWM 27 was calculated by using seven control marks, six 
control marks (excluding B 21 or J 697) and five control marks (excluding B 21 
and J 697). Comparisons of H in GWM 27 are given in Table 33. 


Table 33. COMPARISONS OF H USING FIVE-PARAMETER 
MODEL AT GWM 27 










GWM 27 Hops - Hıeverums (Ст) 
6 (excluding B21) -9.77 -8.07 
6 (excluding J 697) -4,87 -3.65 


b. Heaps on the NPS Campus 
Hops in GH7 and GH8 were calculated by using seven control marks. 
Comparisons of H in GH7 and GH8 are given in Table 34. 
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Table 34. COMPARISONS OF H USING FIVE-PARAMETER 
MODEL AT GH7 AND GH8 


Hon Fomine Cm) 
Check mark 


GH7 0.52 0.62 
GH8 -0.20 0.00 
2. Hops from Six-Parameter Geoid Model 
a. Hops іп the off Campus Area 


Hops in GWM 27 were calculated by using seven control marks, and 


six control marks (excluding B 21 or J 697). Comparisons of H at GWM 27 are 
given in Table 35. 














Table 35. COMPARISONS OF H USING SIX-PARAMETER 
MODEL AT GWM 27 


GWM 27 Hops - Hıeverumg (CM) 
# of control marks 
-12.40 -11.13 


6 (excluding J 697 -3.49 -1.01 
b. Hops on the NPS Campus 


Hops Of GH7 and GH8 were calculated by using seven control marks. 
Comparisons of H at GH7 and GH§8 are given in Table 36. 







Table 36. COMPARISONS OF H USING SIX-PARAMETER 
MODEL AT GH7 AND GH8 


Hos Tiwana Єл 
TEM: 


1.22 1022 
022 0.25 
Since the accuracy of predicted geoid height at check marks (4 0.2 to 2 
cm) is less than the accuracy of GPS observations (+ 2 to 4 cm), it may be 
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concluded that P = I is a satisfactory covariance matrix and further computation 
is unnecessary. 


B. ACCURACY 

The best accuracy of the five-parameter geoid model is € 2 cm in the off 
campus area and is + 2 to 10 mm for the NPS campus. The best accuracy of six- 
parameter geoid model is + 1 cm in the off campus area and is + 2 to 10 mm on 
the NPS campus. 


C. DISCUSSION 
The errors in the geoid model include errors associated with GPS and with 
levelling. 
1. Errors Associated With GPS 
a. Atmospheric Delays 

The time delay of RF signals passing through the ionosphere is due to 
a reduction in speed and the bending of the ray, both effects being due to 
refraction. The overall delay in the signal is nearly inversely proportional to the 
square of the frequency. The ionospheric delay may be calculated by using two- 
frequency receivers (C/A code, P code), and corrected for with a satisfactory 
degree of accuracy. Tropospheric delays are independent of frequency. They 
are relatively small and can be modelled using the elevation angle of the satellite. 
The Trimble 4000SX receiver only observes the C/A code . The Trim640 
program contains no ionospheric modelling. The modified Hopfield model [Fell, 
1975] is used for the tropospheric delay. The combined effects of unmodelled 
ionospheric and tropospheric errors are estimated to result in a satellite-to-user 
range error of from 2.4 to 5.2 m [Milliken, 1980]. 

b. Selection of Appropriate Control Marks 

Only eight permanent bench marks were recovered in the off campus 
study area. D 697 was obstructed by a tree. F 813 was set near a steel bridge. 
Although offset levelling was performed, the ellipsoid height could have been 
affected by offsetting. The shape of the bench marks were spread out roughly in 
a rhomboid shape. B 21 and J 697 are on the ends of the rhomboid. The 
distance is about 33 km from B 21 to J 697. The elevation changes form 5.787 to 
85.061 m. If more bench marks could be recovered and be connected to these 
marks, the accuracy of the seven control marks could be improved. 
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c. Error of Antenna Height Measurements 
The slant heights of GPS antenna were measured by using a Steel tape 
before and after data collection. Thus, the heights of antenna were calculated. 
The systematic error of the tape, reading error, calculation error and a slack tape 
doing measurement will affect the accuracy of the ellipsoid height difference. 
The error of antenna height measurement is + 2 to 10 mm. 
d. Selecting observation period with optimal satellite geometry 
Because there are only seven GPS satellites in operation at present, 
the observing window is only four to five hours a day for four or five satellites. 
The observation schedule should be planed early and considered the travel time, 
and time to reoccupy the marks. In this study most of the GPS measurements 
were made at night. PDOP were not always optimum. 
e. Error in the Computing of Trim640 Program 
There is no standard specification of the GPS software at present. 
Round off error in the data processing, the data record frequency, and the 
precision of the ephemeris could affect the components of the baseline. 
Specifications of the GPS software certainly will be developed in the future. 
2. Errors in Levelling 
The error in levelling would cause the errors of the observation 
equations when determining the local geoid model. The maximum allowable 
closing error between the forward and backward running of a section is 3.0 mm 
Vk for first-order class I levelling, and 9.0 mm Vk for third-order class I 
levelling, where k is the length of the section measured in km [Bodnar, 1975]. 
The average distance from K 152 to the permanent bench marks is about 13 km 
in the off campus study area. Thus, for first-order class I levelling, the 
maximum allowable closing error would be about 11 mm in the off campus area. 
The first-order bench marks were set around 1933. The elevations of the marks 
may have changed caused by subsidences or earthquakes. Since the time did not 
permit rerunning the levelling for the off campus area, published elevations were 
assumed correct. Errors in levelling could contribute to errors in the geoid 
model. 
The distance of the interlocking levelling circuit on the NPS campus is 
about 0.8 km. For third-order class I levelling, the maximum allowable closing 
error is about 7 mm on the NPS campus. The 5 mm closing error was obtained 
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for that region. The geoid height model thus had a better accuracy on the NPS 
campus. 


63 


VIII. NCLUSIONS AND RECOMMENDATION 


Local geoid models for the Monterey, California, area were determined by 
GPS differential positioning using the method of collocation. The local geoid 
models were based on Rapp's 360 degree x 360 order global geoid models 
determined from gravity measurements. Control data were adjusted by least 
squares to solve for the parameters in the local geoid model. Also studied were 
factors that affected the GPS-measured ellipsoid height differences. These included 
(1) comparing GPS differencing solutions, (2) standard error of GPS observations, 
(3) corrections for surface meteorological values, and (4) observation durations for 
GPS. Conclusions are the following: 

]. The accuracy of a local geoid height, N, is the same, whether the NGS or 
the Trimvec version of Rapp's global geoid model is used, although N,(X,Y,Z) 
calculated using Trimble's Trimvec program gives the best results for points on the 
NPS campus, and N,(X,Y,Z) from the NGS program gives the best results for the 
larger study area. For practical use the orthometric heights can be calculated by 
utilizing the local geoid model. 

2. The areas where the geoid is smoothest lead to the best results. This was 
found by comparing the local geoid model used in the larger Monterey Bay area to 
the local geoid model used on the NPS campus. Also the density of control marks on 
the NPS campus was much higher than for the larger study area. If a higher 
accuracy is required in a large area where the geoid may be undulating, more 
control marks are required. These control marks should be strategically located 
throughout the network in order to determine the geoid slope. 

3. The local geoid model does improve the accuracy of Rapp's global geoid 
model locally. The accuracy for Rapp's model is € 2 m in an area 55 km x 55 km. 
The accuracy for the local geoid model is + 2 cm in an area 30 km x 30 km and € 1 
cm in an area 0.5 km x 0.5 km. 

4. A GPS double difference fixed solution gives the best ellipsoid height 
differences, Ah, for baseline lengths up to 30 km. 

5. The standard error of GPS observations of the ellipsoid height difference, 
Ah, 1s + 2 cm to + 4 cm, so there is no need to perform the collocation with a non- 
unit weight matrix. Also it is not possible to reduce the error of the GPS 
observation below + 2 cm. 
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6. There is no need for local surface meteorological corrections to the GPS 
observations. Default meteorological values (1010 millibars, 20° C and 50% 
relative humidity) are sufficient to meet the required accuracy. 

7. The longer the duration of the GPS measurements, the better the result will 
be. The accuracy of the geoid model for the NPS campus, where GPS 
measurements were made over four- to five-hour periods, is better than the 
accuracy of the geoid height model in the larger study area, where measurements 
were made for periods of only ninety minutes. 

8. A five-parameter geoid model gives the best solution for the NPS campus, 
while a six-parameter geoid model gives the best solution for the larger study area. 
The five-parameter geoid model which does not include a AZ term should be used 
with caution where AZ greatly exceeds 5 m. 

Recommendations are as follows: 

1. The GPS antenna should be set at least 1 m above the ground to reduce the 
multiple effects from the ground. The height of the antenna should be carefully 
measured before and after the data collection. 

2. The configuration of the control marks should be carefully planned. Sites 
should be selected to optimize connections to stations with known precise 
orthometric height differences, minimizing the effects of obstructions. Control 
marks in a large area should be numerous enough to meet the required accuracy. 

3. To minimize computation errors least squares adjustments should be made 
using double precision processing when solving the geoid model. For large 
amounts of GPS data a mainframe computer should be used rather than a 
microcomputer to reduce processing time. 
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APPENDIX A. FORTRAN PROGRAM GPSCON 


ok ok ok ok ok ok ok ale ale ale ale ale ale ale ale ade ade ok ale ok ale ok ok ok ok ok ok ok ok ok ale ale ale ale ade ale ade ade ok ok ok ok ok ok ok ok ok ale ale ade ale ale ale ade ade ade ade ade ale ok ae ale ale ale ale ale al ok oes 


WRITTEN BY MA, WEI-MING 05/08/88 
THIS PROGRAM READS DATA FROM GEOID DATA FILE AND CALLS 
THE SUBROUTINE CONTOR TO PLOT GEOID HEIGHTS. 
THE VARIABLES USED ARE: 
GEOHT :  GEOID HEIGHTS (m) 
NX : NUMBER OF POINTS IN THE X-DIRECTION 
NY : NUMBER OF POINTS IN THE Y-DIRECTION 


X 3X eee He * + 
жж SB EF + 


Ae dde ade ale ale ole al ale ale ade ade ade ade ade ale ale ade ale ale ale ale ale ale ale ale ale ale ale ale ale ale ade ade ade ade ae ale ade ade ade ade ade ade ade ade ade ade ale ade ale ale ale ale al ale ale ale ale ade ad ade ale ade ale ade ale al ale ale 


REAL*4 GEOHT(22,15) 

DATA  NX,NY/22,15/ 
C 
C READ GEOID HEIGHTS FROM GEOID1 DATA FILE 
C 


CALL EXCMS(‘FILEDEF 8 DISK GEOID1 DATA A1’) 
DO 200J = 1, NY 
DO 100I = NX, 1, -1 
READ(8,*) GEOHT(LJ) 
100 CONTINUE 
200 CONTINUE 
G 
C PLOT THE GEOID HEIGHT BY CALLING SUBROUTINE CONTOR 
G 
CALL CONTOR(GEOHT,NX,NY) 
STOP 
END 
SUBROUNTINE CONTOR(A,NX,NY) 


ak ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ale ale ale ok ok KK RI KK KK ok ake ake ale ale ale ade ade ade ae ad ade ade ale 


WRITTEN BY MA, WEI-MING 05/08/88 
THIS SUBROUTINE CONTOURS AN NX BY NY ARRAY OF 
REGULARLY SPACED POINTS. 
NOTE: THIS ARRAY MUST BE REAL*4. 
THE VARIABLES USED ARE: 
A : SINGLE PRECISION NX BY NY ARRAY OF 
REGULARLYSPACED POINTS. 
NX : NUMBER OF POINTS IN THE X-DIRECTION 
NY : NUMBER OF POINTS IN THE Y-DIRECTION 
ZINC : CONTOUR INTERVAL 


3€ Y Y Y Ж Ж X EE HE ж 
E A A A A LX A He He He + 


ak ak ak ok ok ok ok ok ole ad ok al ade ade ale ade a oke ok ok ok ok ok cc K 


DIMENSION A(NX,NY) 
COMMON WORK(5000) 


С 
C SET PARAMETERS FOR AXES 
С 
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XORIG = -55.0 
XSTP = 2.0 
XMAX = -34.0 
YORIG = 35.0 
YSTP = 2.0 

Y MAX = 49.0 


€ 
C SET CONTOUR LEVEL 
€ 


ZINC = 0.1 
CALL COMPRS 
CALL SETCLRCCYAN)) 


E 
C SET PAGE AND PLOT SIZES, SET UP AXES FOR PLOT 
С 


CALL PAGE(8.0,6.0) 
CALL BCOMON(5000) 
CALL AREA2D(7.0,5.0) 


С 
C LABEL AXES 
С 


CALL XNAME(LONGITUDE 121 DEGREE WEST (MINUTES)$',100) 
CALL YNAME(LATITUDE 36 DEGREE NORTH (MINUTES)$',100) 
CALL GRAF(XORIG,XSTP,XMAX, YORIG, YSTP, YMAX) 
CALL FRAME 

С 

C MAKE CONTOURS AND DRAW 

С 


CALL SETCLR('RED) 

CALL CONMIN(3.0) 

CALL CONANG(60.) 

CALL CONLIN(0, MYCON',LABELS',1,10) 
CALL CONMAK(A,NX,NY,ZINC) 

CALL CONTUR(,LABELS',DRAW^) 


С 
C END PLOT 

CALL ENDPL(0) 

CALL DONEPL 

RETURN 

END 
C SUBROUTINE MYCON(RARAY,IARAY) 
A de de de dea al al ale ade ade ad ole k ale ol de ad ale ale ale k k e ale ale oie oe oe ale ale ode ode ae a ade ode ale le ae a ode ode je ae ae ade ale ode ade e ae ae ale ale ade ade ade ade ale ale ale ade ok ade K K K 
* ж 
* THIS SUBROUTINE MAKES NEGATIVE CONTOURS DASHED AND * 
* THE ZERO LINE HEAVYIER : 
* ж 
A od od de K k k k k k ade ae ode ade k k ode ode a al de a ode ode oe ade ode ode ale ale ade a ade ade a ade ade ale ale de a ale ode a ade ode ade a ade ode ale ade ode ade ae ade ade ade ae ade ade k k kK k K ak 

DIMENSION RARAY(2),IARAY (1) 

CALL RESET(DASH) 

IF (RARAY(1).GE. 0.) GO TO 10 

CALL DASH 

10 RARAYQ)=1. 
IARAY(1) =1 


67 


IF (RARAY(1).EQ.0.) IARAY(1) =2 
RETURN 
END 
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356 
357 
358 
359 
380 


APPENDIX B. BASIC PROGRAM LOBS.BASIC 


REM LEAST SQUARES BY OBSERVATION 
DIM L(N%,1) 

DIM A(N%,X%), P(X%,X%), EX(X%,1) 
DIM CX(X%,1) 

DIM AT(X%,N%), R(50), АТР(Х%,Х%) 
DIM ATPA(X%,X%), ATPL(X%,1), AUX%,X%) 
PRINT " INPUT # OF OBSERVATIONS" 
NPUT NO 

PRINT " INPUT # OF PARAMETERS" 
INPUT N1 

ERASE L 

ERASE A 

ERASE AT 

ERASE P 

ERASE ATP 

ERASE ATPL 

ERASE AI 

ERASE CX 

ERASE EX 

ERASE ATPA 

DIM L(NO,1), A(NO,N1), AT(N1,NO), P(NO,NO), ATP(N1,NO), 


ATPL(N1,1),AIC(N1,N1) 


390 
400 
1180 


DIM CX(N1,1), EX(NO,1), ATPA(N1,N1) 
ITERO = 0 
PRINT " INPUT COEFFICIENTS OF OBSERVATION EQUATIONS AND 


WEIGHTS" 


1190 
1200 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
1300 
1390 
1410 
1420 
1430 
1450 
1460 
1470 
1480 
1490 
1500 
1510 
1520 
1530 
1540 


FOR K = 1 TO NO 

FOR J-1TONI 

PRINT "COEFFICIENT ", К, Ј 

INPUT A(K,J) 

PRINT A(K,J) 

NEXT J 

PRINT " INPUT OBSERVED VALUE ", K 
INPUT L(K,1) 

PRINT " INPUT WEIGHT OF OBSERVED VALUE ", K 
INPUT P(K,K) 

NEXT K 

DOF =0 

REM LEAST SQUARES 

М = МО 

І. = № 

РВІМТ "М =", М, "№ =", №, "= "1. 
FORI=1TOM 

FOR J=1TOL 

АТ(Ы) - А(1,7 

NEXT J 

NEXT I 

REM ATP = AT *P 

N=NO 

FOR I=1 TOL 
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1550 
1560 
1570 
1580 
1590 
1600 
1610 
1620 
1630 
1640 
1650 


1680 
1690 
1710 
1720 
1730 
1740 
1750 
1760 
1770 
1780 
1790 
1800 
1810 
1820 
1830 
1840 
1970 
1980 
1990 
2000 
2020 
2030 
2040 
2050 
2055 
2060 
2070 
2080 
2090 
2100 
2110 
2120 
2130 
2140 
2150 
2160 
2170 
2180 
2190 
2200 
2210 


FOR J=1 TON 

ATP(LJ) =0 
FORK=1 TOM 
ATP(LJ) = ATP(,J) + AT) * P(K,J) 

NEXTK 

NEXT J 

NEXT I 

REM ATPA - ATP* A 

N-NI 

FOR I=1 TOL 

FOR J=1 TON 

ATPA(1) - 0 

ATPA(LJ) = 0 

FOR K=1TO M 

ATPA(L, J) = ATPAA,J) + ATP(LK) * A(K,J) 

NEXTK 

NEXT J 

NEXT] 

REM ATPL- ATP*L 

N=1 

FORI=1TOL 

FORJ=1TON 
ATPL(LJ) = 0 
FOR K=1 TOM 
ATPL(I,J) = ATPL(I,J) + ATP(I,K) * L(K,J) 

NEXT K 

NEXT J 

NEXT I 

REM AI - INV(ATPA) 

I-NI 

М = № 

N=I-1 

МІ= М -1 

ЕОВ Ј = 1 TOI 

FOR K =1 TOI 

AIQ,K) = ATPA(J,K) 

NEXT K 


NEXT J 
FOR K -1 TOI 
FOR J = 1 TO MI 

R(J) = AI(1,J+1) / AI(1,1) 
NEXT J 

R(M) = 1! / AI(1,1) 
FOR L=1 TON 
FOR J = 1 TOMI 


AI(I,J) = AI(L+1,J+1) - AI(L+1,1) * R(J) 
NEXT J 
AI(L,M) = -AI(L+1,1) * R(M) 
NEXT L 
FOR J = 1 TO M 
AI(I,J) = R(J) 
NEXT I 
NEXT K 
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2220 
2230 
2240 
2250 
2260 
2270 
2280 
2300 
2310 
2320 
2325 
2330 
2340 
2480 
2490 
2500 
2510 
2520 
2530 
2540 
2560 
2570 
2580 
2590 


2610 
2620 
2630 
2640 
2650 
2660 
2670 
2680 
2690 
2695 
2700 
2710 
2720 


REM CX = AI * ATPL 


L=NI1 

N =1 

MzNI 

FOR I=1 TOL 
FOR J =1 TO N 
CX(I,J) = O 
FOR K=1TOM 


CX(1J) = CX(1,J) + AI(LK) * ATPL (K,J) 
NEXT K 
PRINT "PARAMETER # ", I, CX(LJ) 


NEXT J 

NEXT I 

PRINT " RESIDUAL" 
L=NO 

N=1 

MzNI 
REM EX = A * CX 
FOR I=1 TO L 
FOR J =1 TON 
EX(1,J) = 0 

FOR K=1 TOM 


EX(I,J) = EX(1,J) + ACK) * CX(K,J) 
NEXT K 

NEXT J 

NEXT I 

SD =0 

FOR K =1 TONO 

V = EX(K,1) - L(K,1) 

PRINT K, V 

SD=SD+V*V 

NEXT K 

SD = SQR(SD / (NO - N1 + DOF) 
PRINT SD 

PRINT " STD.DEV", SD 

ITER = ІТЕВ + 1 

IF ITER <3 GOTO 350 

END 
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APPENDIX C. FORTRAN PROGRAM GPSDIS 


KK K K K K ok oko ok ok ok ok ok ok ok ok oko ok ale ale ale ak ade al ale ale ale ae ae ade ade ale ale ale ade ak ak ale ale ale ak ade ad oko ale ale ae ade ade ak ad ade ade ade ale ade ad ade ale ale ak ade ok al ak al ak ak 


* WRITTEN BY WEI-MING MA 05/14/88 ж 
* — THIS PROGRAM READS DATA FROM TERMINAL AND COMPUTES ж 
*  THEORTHOMETRIC HEIGHTS FOR THE CHECK MARKS THEN ж 
* COMPUTES THE DIFFERENCE. * 
* THE VARIABLES USED ARE: ж 
ж RX,RY,RZ : THE COORDINATES OF THE REFERENCE MARK (m) * 
* CX,CY,CZ : THE COORDINATES OF THE CHECK MARK (m) * 
* DX,DY,DZ : THE COORDINATES DIFFERENCE (m) * 
ж Н : THE ORTHOMETRIC HEIGHT OF CHECK MARK (m) * 
ж LH  : THELEVELED ORTHOMETRIC HEIGHT (m) ж 
ж DH : THE ELLIPSOID HEIGHT DIFFERENCE (m) * 
* GN  : THE GLOBAL GEOID HEIGHT OF CHECK MARK (m) * 
* ІН : THE LEVELED ORTHOMETRIC HEIGHT (m) * 
ж DIF : THE DIFFERENCE BETWEEN H AND LH (cm) ж 
ж A,B,C,D : THE COEFFICIENTS OF 5-PARAMETER GEOID MODEL * 
ж A1,B1,C1 : THE COEFFICIENTS OF 6-PARAMETER GEOID MODEL * 
ж D1,E1 ж 
ж IPAJPA1 : THE NUMBER OF PARAMETERS * 
ж IW : THE OUTPUT UNIT ж 
ж ж 


OO I I IGG ЖЖЖЖ ЖЖЖ KK aK 
REAL*8 RX, RY, RZ, CX, CY, CZ, DX, DY, DZ, LH, GN, 
REAL*8 A, B, C, D, A1, B1, C1, D1, El, H, DH, DIF, 
INTEGER IPA, IPA1, IW 
CHARACTER*1 RESPN1, RESPN2, RESPN3 
IW-6 
IPA =0 
IPA1 20 


С 
C READ DATA FROM TERMINAL 
G 
10 PRINT *, ENTER # OF PARAMETERS' 
READ *, IPA 
IF (IPA .EQ. 5 .OR. IPA .EQ. 6) THEN 
PRINT *, 'ENTER THE X COORDINATE IN WGS 84 OF THE ', 
& ‘REFERENCE MARK' 
READ *, RX 
PRINT *, 'ENTER THE Y COORDINATE IN WGS 84 OF THE ', 
& "REFERENCE MARK' 
READ *, RY 
PRINT *, ENTER THE Z COORDINATE IN WGS 84 OF THE ', 
& REFERENCE MARK' 
READ *, RZ 
ELSE 
STOP 
END IF 
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IF (RESPN1 .EQ. 'Y') THEN 
PRINT *,'DO YOU WANT TO CHANGE THE CHECH MARK ("Y" OR ', 
& rs "М et ) t 
READ *, RESPN3 
IF (RESPN3 .EQ. "Y') GO TO 30 


GO TO 40 
END IF 
20 PRINT *,'DO YOU WANT TO CHANGE THE REFERENCE MARKS ("Y" ', 
& 'OR UNT) 


READ *, RESPNI 
IF (RESPNI .EQ. 'Y') GO TO 10 
30 PRINT *,ENTER THE X COORDINATE IN WGS 84 OF THE CHECK. 
& 'MARK' 
READ *, CX 
PRINT *,,ENTER THE Y COORDINATE IN WGS 84 OF THE CHECK MARK’ 
READ *, CY 
PRINT *,,ENTER THE Z COORDINATE IN WGS 84 OF THE CHECK MARK' 
READ *, CZ 
PRINT *,'ENTER THE ELLIPSOID DIFFERENCE, DH’ 
READ *, DH 
PRINT *,ENTER THE GLOBAL GEOID HEIGHT OF THE CHECK MARK; 
& " GN' 
READ *, GN 
PRINT *,'ENTER THE LEVELED OTHOMETRIC HEIGHT OF THE CHECK ' 
& , MARK' 
READ *, LH 
DX = CX - RX 
DY =CY -RY 
DZ = CZ- RZ | 
IF (IPA .EQ. IPA1 .AND. RESPNI .EQ. 'N' .AND. IPA .EQ. 5) GO TO 50 
IF (IPA .EQ. IPA1 .AND. RESPNI .EQ. 'N' .AND. IPA .EQ. 6) GO TO 60 
40  IF(IPA .EQ. 5) THEN 
PRINT *,ENTER THE ELLIPSOID HEIGHT H0' 
READ *, HO 
PRINT *,ENTER THE COEFFICIENT A' 
READ *, A 
PRINT *,ENTER THE COEFFICIENT B' 
READ *, B 
PRINT */ENTER THE COEFFICIENT C' 
READ *,C 
PRINT *,ENTER THE COEFFICIENT D' 
READ *, D 


C 
C THE 5-PARAMETER GEOID MODEL 
e 


50 H = -GN + DH + HO + A*DY + B*¥DX**2 + C*¥DY**2 + D*¥DX*DY 
ELSE 

PRINT *,, ENTER THE ELLIPSOID HEIGHT НО' 

READ *, HO 

PRINT *,'ENTER THE COEFFICIENT A" 

READ *, Al 

PRINT *,'ENTER THE COEFFICIENT B" 

READ *, B1 
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PRINT *, ENTER THE COEFFICIENT C" 
READ *, Cl 
PRINT *,'ENTER THE COEFFICIENT D" 
READ *, D1 
PRINT *,'ENTER THE COEFFICIENT E" 
READ *, El 


С 
C THE 6-PARAMETER GEOID MODEL 
С 


60 H = -GN + DH + HO + Al*DX + B1*DZ + C1*DX**2 + DI*DY**2 
& + El*DX*DY 
END IF 
С 
C COMPUTE THE DIFFERENCE AND PRINT THE RESULT 
С 


DIF = (H - LH) * 100. 
WRITE(IW, 1) IPA 
WRITE(IW, 2) LH, H, DIF 
1 FORMAT (5X, I1,' PARAMETERS) 


2 FORMAT (/5X,H-LEVELING -',F10.3,' M' 
& /SX,H-GPS =', F10.3, ' M' 
& /5X,DIFFERENCE =’, F10.3, ' СМУ) 
ІРА1 -ІРА 
Н - 0. 
DIF = 0. | 


PRINT *, MORE COMPUTATION ("Y" OR "N")' 
READ *, RESPN2 

IF (RESPN2 .EQ. 'Y') GO TO 20 

STOP 

END 
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APPENDIX D. FORTRAN PROGRAM DISTCO 


ЖЖЖЖЖЖЖЖЖЖЖЖҰЖЖЖЖЖЖЖЖҰЖЖЖЖ ЖЖ ЖЖ ЖЖЖ ЖҰ ЖЖЖЖЖЖЖЖ ЖЖЖ ЖЖЖ ЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖ 


ж — WRITTEN BY MA, WEI MING 05/11/88 Š 
ж THIS PROGRAM COMPUTES THE DISTANCE BETWEEN THE BENCH * 
* МАКК THEN CALLS THE SUBROUTINE C TO CALCULATE THE ы 
* COVARIANCE BETWEEN THE BENCH MARKS. 5 
ж THE VARIABLES USED ARE: * 
E X,Y,Z : THE COORDINATES OF X,Y,Z, IN WGS 84 (m) К 
* DIS : THE DISTANCES BETWEEN THE MARKS (m) Ы 
* CQ :  THECOVARIANCES BETWEEN THE BENCH MARKS * 
* * 
ж 


ЖЖ ЖЖ ЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖЖ ЖЖЖ ЖЖЕЖЖЖЖКЖЖЖЖЖЖЖЖЖЖЖЖЖ ЖЖЖ ЖЖЕЖ ae aaa aaa aaa a 
REAL*8 X(10), Y(10), Z(10) 
REAL*8 DIS(10,10), CQ(10,10) 
DATA  DIS/100*0/, CQ/100*0./ 
IW =9 


N=7 
С 
C READ X, Y, Z FROM THE DATA FILE 
€ 


CALL EXCMS(FILEDEF 8 DISK MBXYZ DATA Al) 
DO 100I2 L, N 
READ(8,*) XD, YO, ZU) 
100 CONTINUE 


E 
C COMPUTE THE DISTANCE AND PRINT THE RESULTS 
e 


DO 300 1=1,N 
WRITE(W,1) I 
1 FORMAT(/1X,FROM BENCH MARK 2,12, TO', DISTANCE (m)/ 
& 1Х, 53(-')) 
ро 200] = 1, № 
IF (J EQ. 1) GOTO 10 
DIS(I,J) = SQRT((X(J) - X(1)) ** 2 + (Y(J) - Y(1)) ** 2 


& +(Z(J)-Z(1))**2) 
10 WRITE(IW,2) J, DIS(I,J) 
р FORMAT(11X, 12, 14X, G20.8) 


200 CONTINUE 
300 CONTINUE 


С 
C CALL SUBROUTINE C TO CALCULATE THE CQ AND PRINT THE RESULTS 
e 


WRITE(IW,3) 
3 FORMAT(/3X, 'CQ =') 
CALL C(DIS,N,CQ) 
DO 400J =1,N 
WRITE(IW,4) (CQ(LJ), I = 1, №) 
4 FORMAT(3X, 7(G16.7, 2X) 
400 CONTINUE 
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STOP 
END 
SUBROUNTINE C(DIS,N,CQ) 


CE E E E EE EE E E E EEEE EEEE EEEE EE EE E E EE EEEE EEE EEE EEEE E E E E E E E E E E E E E E E E E E E EE E E E 
ж 


Ж 
* — WRITTEN BY MA, WEI MING 05/11/88 ж 
* THIS SUBROUTINE COMPUTES THE COVARIANCES OF THE RANDOM * 
* ERROR DISTRIBUTION FUNCTION: 

*  C(R)=A+B*SIN(D*DIS) 

* THE VARIABLES USED ARE: 

ж А : STANDARD DEVIATION OF THE OBSERVATIONS 
ж В : CONSTANT 

ж С : CONSTANT (RADIANS/METERS) 

ж DIS  : DISTANCE (METERS) 

ж 

Ж 


+ << + <= © + 


ЖЖЖЖЖЖЖ ЖЖЖЖ 


REAL*8 DIS(10,10), CQ(10,10), A, B, D 

А = 0.137153 

В = 0.147951 

D = 2.85958 

DO 200 J=1,N 

DO 1001=1,N 
CQ(LJ) 2 A - B * SIN(D * DIS(IJ) / 33300.) 

100 CONTINUE 
200 CONTINUE 

RETURN 

END 
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